Nonlocal density functionals and the linear response of the homogeneous electron gas 
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The known and usable truly nonlocal functionals for exchange-correlation energy of the inhomo- 
geneous electron gas are the ADA (average density approximation) and the WDA (weighted density 
approximation). ADA, by design, yields the correct linear response function of the uniform electron 
gas. WDA is constructed so that it is exact in the limit of one-electron systems. We derive an 
expression for the linear response of the uniform gas in the WDA, and calculate it for several flavors 
of WDA. We then compare the results with the Monte-Carlo data on the exchange-correlation local 
field correction, and identify the weak points of conventional WDA in the homogeneous limit. We 
suggest how the WDA can be modified to improve the response function. The resulting approxi- 
mation is a good one in both opposite limits, and should be useful for practical nonlocal density 
functional calculations. 



Calculations based on the Kohn-Sham formulation of 
density functional theory^ have become a prominent tool 
in condensed matter physics. Current work is dominated 
by local density approximation (LDA) studies, in which 
the exchange correlation functional is a local function of 
the density. However, as the number and accuracy of cal- 
culations has increased, so has the number of well doc- 
umented cases where the LDA is inadequate and with 
this interest in beyond LDA approaches, e.g., the gen- 
eralized gradient approximation (GGA), which depends 
locally on both the density and its gradient. 

Modern GGA functionals do improve upon LDA re- 
sults for a wide range of problems. However, several 
studies have pointed out deficiencies in GGA function- 
als, e.g. difficulties in describing ferroelectric materials, 
and cases of overcorrection of LDA errors particularly in 
materials containing heavy atoms. Both the LDA and 
GGA fail to provide a correct description of the static 
short-range linear response of the homogeneous electron 
gas. All this leads to the question of the extent to which 
truly non-local functionals are practical and able to cor- 
rect the deficiencies of LDA and GGA methods. 

The first efforts at developing practical non-local func- 
tionals date from the 1970's when the weighted den- 
sity approximation (WDA) and average density approx- 
imation (ADA) were proposed. However, over most of 
the intervening period the field has been relatively dor- 
mant, in part because of the success of the simpler LDA 
and GGA schemes and in part because it was widely 
thought that such schemes could not be implemented in 
a computationally tractable fashion. However, at least 
for the WDA, computationally efficient algorithms are 
now known^"^ and benchmark calculations have been re- 
ported. In the cases that have been studied ground state 
properties are generally improved over the LDA^ . 

Both methods were proposed in 1978-1979 by Alonso 
et al and Gunnarsson et af'"^. Both exploit the general 
expression for E^c, 



n(r)n(r') 



G(r, r'){n{r)}drdr' , 



(1) 



where the function G(r, r') is also a functional of 
the total electronic density n{r). A rigorous expres- 
sion for G can be derived^ in terms of coupling con- 
stant averaged pair correlation function: G(r, r') — 
[g(r, r'; A){n(r)} — l]dA. For the uniform gas this func- 
tion, Go(|r — r'l, n), is known with high accuracy^, but 
for an arbitrary system there is no practical way to 
use this formula. The LDA instead of Eq. (1) uses 
(eV2)/(irrfr'n2(r)Go[|r-r'|,7i(r)]/|r-r'|, so that E^^ 
becomes E^^^ = J n{r)€xc[n{r)]dr , e^c being the 
density of exchange-correlation energy of the uniform 
gas. The LDA is incorrect in the two important lim- 
its: the fully localized, i.e., a one electron system, and 
the fully delocalized limit, i.e., homogeneous electron 
gas. In the former case the LDA gives a spurious self- 
interaction with energy (e^/2) J (ir(ir'n(r)7i(r')/|r — r'|-|- 
/ n{Y)€xc[n{v)]dY , which is widely thought® to be a key 
problem with the LDA. In the homogeneous limit, the 
LDA gives the correct exchange-correlation energy, but 
the changes of this energy upon small perturbations are 
not properly described; the second variation of E^c with 
density, i.e., the exchange-correlation part of the dielec- 
tric response, K^d^ — r') = 5"^ E^d 5'n{v)8n{Y'), is a delta 
function, which is incorrect. The Fourier transform of 
Kxcif) in LDA is independent of the wave vector. Since 
LDA is exact for the uniform gas, K^^^ corresponds to 
the correct K^c at g = 0. GGAs also give correct behav- 
ior at g = 0, but become even worse than the LDA at 
high g's. 

The two nonlocal expressions for E^c, WDA and ADA, 
aimed at correcting one or the other of these two limits. 
The former uses the general expression (1), but instead 
of the actual function G uses a model function, defined 
so that the one electron limit is honored. This begins 
by choosing a generic expression for G, which depends 
on one parameter n, to be defined later. In the original 
papers it was suggested that G(r, r',n) = Gh{'r,r',h) = 
J^[g{r,r'; X,n) — l]dX., where g is the pair correlation 
function of the homogeneous electron gas. Later it was 
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realized-'^^ that other choices of G may be better than Gh- 
In the WDA n is a function of r, but differs from n(r), 
and is chosen so that / G[r,r' ,n{r)]n{r')dr' = —1. This 
assures that for a one electron system E^c cancels the 
self-interaction exactly. 

In the ADA n(r') in Eq. (1) is substituted by n(r), 
which results in E^^^^ = / n(r)e2;c[r7,(r)]c?r. Then n(r) is 
defined as n(r) = Jii;[|r — r'|, n(r)]n(r')(ir', and the uni- 
versal function w is chosen so that 5'^ E^^"^ / 5n{Y)5n{v') 
gives the correct K^c for the uniform gas. Contrary to 
the WDA, the ADA is not self-interaction free in one 
electron systems. 

From the beginning there was substantial interest in 
the behavior of WDA in the dclocalizcd limit^. Williams 
and von Barth^^ suggested that the WDA should give 
substantial improvement over the LDA in this limit, but 
till now no systematic study has been reported. If this 
conjecture is true, the WDA has a great advantage over 
any other known approximation to the DFT in the sense 
that it accurately reproduces two key physical limits. 
Furthermore, even if it is not entirely correct, the next 
question is, whether or not an approximation based on 
the WDA exists that does provide proper limiting behav- 
ior. In this paper we derive an expression for Kxc in the 
WDA, calculate Kxc for popular flavors of WDA, and 
discuss construction of an improved WDA method. 

We start by deriving a closed expression for K^c in the 
WDA for an arbitrary G. First some notation: we denote 
the product {e^ /r)G{r) as W{r), use atomic units where 
e = 1, h ~ 1, and use primes for the derivative with re- 
spect to the density argument, e.g. G' = dG/dn. We also 
introduce two functions, reflecting implicit dependence of 
the weighted density n on variations of the real density: 

d{r'-r) = Snir')/Sn{r) (2) 
S-ii(r') Sd{r'-r) 



/(r'-r,r'-r") 



6n{r)dn{r") (5n(r") ' 



(3) 



Using the WDA expression for the exchange-correlation 
energy, 

E,, = (1/2) J n{v)n{r')W[\v - r'\,n{v)]drdv' , (4) 

we can express Kxc in terms of functions d and /, and we 
can find these functions using normalization condition 



dr'n{r')g[\r — r'|, n(r)] 
and the LDA limit condition 

dr'W[\r-r'\,n\=2exc/n. 



1, 



(5) 



(6) 



We proceed in reciprocal space, which corresponds to us- 
ing density perturbation of the form i5n(r) =nqe^'^^ . Let 

Wq, Gq, dq and fp q will be the Fourier transforms of the 
corresponding functions. Then the above condition can 
be written as 



Go = -l/n, Wo = 2exc/n. (7) 

For dq one finds dq = —nGq. For fp,q we need only di- 
agonal elements, fq,-q = 2nGq{nG'q + Gq). In terms of d 
and /, Kxc is 

Kxc{q) = Wq+nodqW:j+nodqW;,+^{dlw;,'Wofq,-qW^). 
resulting in 

Kxc{q) =Wq- n^Gq{W'q + VF^) + n'{n'Glw;,)'/2. (8) 
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FIG. 1. Exchange-correlation local field factor in the 
WDA of Ref. 10 (Gunnarsson-Jones), Ref. 12 (Gritsenko 
et al), and derived from the homogeneous electron gas 
pair correlation function (Perdew-Wang), as compared 
with the Monte Carlo results (Monte Carlo) and the in- 
terpolating formula thereof (MC interpolation), as given 
in Ref. 13. Densities, from top to bottom, correspond to 
Ts = 1,2,5. 
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The original formulation of the WDA used the cor- 
responding homogeneous electron gas function for G. 
Since then, three forms of G have been used in the 
calculations, all of which result in improvement over 
LDA (in the admittedly limited number of tests per- 
formed to date). These are: the function G de- 
rived for the uniform gas by Perdew and Wang^, 
the Gunnarsson- Jones function G'^'^ {r) = Ci(n){l — 

exp[— ^ C2(n) ) ]}' ^^'^ Gritsenko et al^^ function 

qGrba^^^ = Ci(n)exp[- (^-^)'']}, k = 1.5 (note that 

the uniform gas function^ is approximately given by the 

same expression with k = 2). We tested these func- 
tions for the densities = 1, 2, and 5 and obtained 
modest agreement with the Monte Carlo results^^ (Cf. 

Fig.l, where we plotted calculated exchange-correlation 

2 

local field factor Ixcil) = 4:;^Kxc{q), and compare it with 
Monte Carlo data^^). By construction, Kj;c{0) is cor- 
rect (and is in fact the LDA value). At 5 > 1.5 — l.Skp 
Kxc falls below its LDA value and continues to decrease 
at large g's. However, a closer look reveals two ma- 
jor disagreements: first, I^^"^{q) is considerably larger 
than the Monte-Carlo data for the wave vectors between 
« O.Sfci? and l.Sfci?. Second, Ixc{q) in WDA tends to a 
constant value equal to lim.q^oo[Wq/ {Air / q^)]. In Monte 
Carlo calculations it is K^dq) itself that has a finite 
limit at q ^ 00, while Ixdq) const ■ q'^ at q ^ 00. 
The latter result was predicted by Holas^'' and is physi- 
cally important: it reflects the fact that E^c is not solely 
an interaction energy, but also includes the exchange- 
correlation contribution to kinetic energy (which in fact 
decays slower with q than the interaction part of Exc)- 



Can one correct these two deficiencies without compro- 
mising the correct one-electron limit of WDA? In fact, it 
was noticed long ago^" that there is no particular rea- 
son to use the homogeneous electron gas pair correla- 
tion fimction for G. Since using Gh in WDA does not 
guarantee any improvement in describing properties of 
the homogeneous gas itself, one may use the freedom in 
G{r) to adjust the WDA so that the calculated local field 
factor (and thus linear response function) is as accurate 
as possible. Inversion of eq.(8) yields G{q) for a given 
Kxc{q)- It does not guarantee, however, that the result 
will be physical. So, as a first step, let us analyze Eq. 
(8). For this purpose, we write Gq = —(p{p/Q)/n, with 
the condition ip{Q) = 1, where Q is some constant (both 
the Gunnarsson-Jones and the Gritsenko et al functions 
are of this form). Then 



d^q- 
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Gqdq = - 
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If we now define Q{n) = —Trexcin), then the second 
condition on ip{x) becomes ip{x)dx = 1. These two 
conditions reduce our freedom to adjust Gq : since the 
characteristic size of ip{x) is of order of 1, the wave vector 
dependence of Gq is defined by the ratio q/Q = —qjixtxc- 
At high density Q — l.SSSfci?. and only close to > 6 
does it approach 2^^^, the number at which real local 
field factor changes its behavior from low-g to the high-g 
limit. A monotonic function (^(x) does not reproduce this 
feature, which explains why existing WDA parametriza- 
tions put the hump in Kxc at too low q. Nonmonotonic 
and explicitly density-dependent functions <^(x) may be 
able to shift the hump to its correct position at q = 2kp. 
It is still an open question whether or not a physically 
sound function can be found with this property. 

However, even if the "2fci?" problem is fixed, another, 
probably even more important problem remains: the 
short wave length behavior of Kxc- Fortunately, this is 
easy to correct. Farid et al^^ tabulated the coefficient 7 
that defines the asymptotic behavior of Kxc{q — > 00) as 

Kxc{q — * 00) = — ■^7(n)|5-- These values can be fit as 

7W-(^j 15 ^^''>~ [l+cx + dx^) 

where a = 0.026319, h = 0.00823859, c = -0.173199, 
d = 0.233081. Let us now modify the function G{r) 

G{r) = Gi{r) + G2{r) = AS{r)/Anr + G2{r), 

Since J Gi{r)r^dr = 0, the normalization condition for 
G2 is the same as for G itself. Since An J Gi {r)rdr = A, 
the LDA limit condition for G2 becomes 

47r j G2{r)rdr = 2ixc{n)/n, 

ixc{n) = exc{n) - An/2. 

Thus 

A = = -(^)^/^^/(^/^) = -3.0856r^/(^^^), 

0.368317 

exc{n) = exc[n) H /(v^) 



Now Gp = G2p, Wp = A + W2p, and W; = A' + W^p, 



Kxc{q) =A + W2q- n^G2q{A' + W^q + Wi^o) + n^n^GlWi^oY = A- nlG2pA' + Kxc, (9) 

where Kxc is calculated from Ixc in exactly the same way as Kxc is calculated from Cxc- The corresponding functional 
for the exchange-correlation energy is 
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E: 



•AWDA 



n(r)n(r') 



W{\t- 



iyG[|r-r'|,n(r)]drdr' + j n{v)^-^^f 



dr. 



(10) 



Here 47rrf/3 = n, and G{r) is normalized to €xc{n) 
instead of exc{n). Since we do not require that G{r) = 

lo [di'^'j ^) ~ i]dX., where g corresponds to the uniform 
gas, but rather consider it to be a flexible function satis- 
fying two normalization conditions, further improvement 
of the method should be possible along the line described 
in the previous paragraph, namely the freedom in chosing 
G{r) can, and should, be used to yield K^c according to 
Eq. (9) close to the linear response of the homogeneous 
electron gas, including correct behavior near q = 2kp. In 
Fig. 2, we show Ixc calculated according to Eq. 10 with 
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the different functional form of G{r). Clearly, the re- 
sults are much better than either the LDA or "con- 
ventional" WDA. Interestingly, when the nearly exact 
Perdew-Wang function, or exponential function with k = 
2, are used, the resulting Ixdl) is close to the analytical 
function derived by Farid et al (arguably the best an- 
alytically derived Ixc{q) available), while an exponential 
function with k = 1.5 is close to the formula of Ref. 13, 
which is a fit to the Monte Carlo data. 

To summarize, wc have calculated the exchange- 
correlation local field function Kxc in the WDA, and 
found that besides the expected improvement over the 
LDA it has two major deficiencies: (1) it does not have 
correct asymptotic behavior at g ^ oo, and (2) the char- 
acteristic feature at g « 2A:f is displaced towards smaller 
q's. The former can be easily corrected by adding a delta- 
function component to G(r), which results in Eq. (10). 
The latter is harder to fix, but there arc still unused de- 
grees of freedom in the formalism which may be used 
to tune the behavior near 2kF- In our opinion, this new 
scheme for WDA calculations (Eq. 10) is currently most 
promising for practical applications. 

This work was supported by ONR. 



FIG. 2. 7a;c(g) as in Fig.l, but for the modified WDA 
of Eq.lO. Also the analytical formula of Farid et al'^^ is 
shown. 



^ Theory of the inhomogeneous electron gcis, Ed. by S. 

Lundqvist and N.H. March, Plenum, 1983. 
^ D.J. Singh, Phys. Rev. B, 48, 14099 (1993), and Ferro- 

electrics, 194, 299 (1997). 
^ M. Sadd and M.P. Teter, Phys. Rev. B, 54 13643 (1996). 

J.P.A. Charlesworth, Phys. Rev. B, 53, 12666 (1996). 

* Ref. 4 is an exception, but note that the LDA results of 
this work differ from most previous calculation. 

^ J.A. Alonso and L.A. Girifalco, Phys. Rev. B, 17, 3735 
(1978). 

O. Gunnarsson, M. Jonson, and B.I. Lundquist, Phys. Rev. 
B, 20, 3136 (1979). 

* J. P. Perdew and Y. Wang, Phys. Rev. B 46, 12947 (1993). 
® See, e.g., S. Goedeckcr and C.J. Umrigar, Phys. Rev. A 55, 

1765 (1997), and references therein. 

O. Gunnarsson and R.O. Jones, Physica Scripta, 21, 394 
(1980). 

" A.R. Williams and U. von Barth, in Ref. 1. 

O.V. Gritsenko, A. Rubio, L.C. Balbas, and J.A. Alonso, 

Chem. Phys. Lett 205, 348 (1993) 
" S. Moroni, D.M. Ceperley, and G. Senatore, Phys. Rev. 

Lett 75, 689 (1995). 
^'^ A. Holas, in StTongly Coupled Plasma Physics, ed. by F.J. 

Rogers and H.E. De Witt (Plenum, NY, 1987) 

B. Farid, V. Heine, G.E. Engel, and I.J. Robertson, Phys. 

Rev. B 48, 11 602 (1993).4 



4 



